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We present a simple density functional theory for the solid phases of systems of particles interacting 
via soft-core potentials. In particular, we apply the theory to particles interacting via repulsive point 
Yukawa and Gaussian pair potentials. We find qualitative agreement with the established phase 
diagrams for these systems. The theory is able to account for the bcc-fcc solid transitions of both 
systems and the re-entrant melting that the Gaussian system exhibits. 

PACS numbers: 



I. INTRODUCTION 

When an ensemble of particles is at a sufficiently high 
density, the fluid will often freeze to form a crystalline 
solid. This freezing is driven by the repulsion that invari- 
ably exists between particles when at close separations. 
The question that naturally arises is: how much repul- 
sion is needed for freezing and what crystal structures 
are formed? The answer is that even particles interact- 
ing via repulsive pair potentials v(r) that are finite for 
all separations r can freeze 

QUO!. The particular crys- 
tal structure that is formed depends upon the form of the 
particle interactions and especially on the 'softness' of the 
decay of v{r) as r — > oo 4]. Furthermore, there can be 
solid-solid transitions, with different crystal structures 
being stable in different portions of the phase diagram. 
In this paper we present a simple density functional the- 
ory (DFT) jjj for determining the location of the melting 
and solid-solid phase boundaries for soft-core particles. 

For the purposes of this paper, we define 'soft-core' 
particles as those with purely repulsive pair potentials 
v(r), for which the integral over all space of v(r) is fi- 
nite, i.e. J drv(r) < oo. Alternatively, one can define 
soft-core potentials as those for which the Fourier trans- 
form v(k) of the pair potential exists. Commonly encoun- 
tered (model) potentials in the theory of liquids such as 
the Coulomb potential, the hard-sphere potential or the 
Lennard-Jones potential |(| do not fall into this cate- 
gory. However, a wide class of fluids can be modelled by 
particles interacting via potentials that do fall into the 
soft-core category. One common example is the Yukawa 
core model (YCM): 

eexp(-Ar) 
v{r) = , (1) 

where A > and e > 0. Such a potential is used to 
model the effective interaction between charged point 
particles, where the Coulomb interaction between the 
particles is screened by a background medium. The ef- 
fects of the screening are incorporated in the parameter 



A [3, |8j . Examples of systems where the particle interac- 
tions can be modelled by a Yukawa pair potential range 
from charged colloidal solutions 0, H E3 to dusty plas- 
mas [Illiallilii. 

Another soft-core potential is the Gaussian core model 
(GCM): 

v(r) = eexp(-AV). (2) 

The freezing behaviour of a fluid composed of particles 
interacting via such a potential aroused much interest 
due to the novel re entrant melting behaviour: for cer- 
tain temperatures, on increasing the fluid density, the 
fluid freezes. However, on further increasing the density, 
the crystal re-melts The high density phase of the 
GCM is the fluid state 0, H El El E E3 . A Gaus- 
sian potential is used to model the effective interaction 
between the centres of mass of polymers, star-polymers 
and dendrimers in solution [lj,|20(. In this case A -1 ~ R, 
the radius of gyration of the polymers. 

The DFT theory for freezing presented here is a simple 
qualitatively accurate theory which, for the cases we have 
tested, gives the correct topology of the phase diagram, 
including the existence of solid-solid phase transitions - 
i.e. the theory incorporates in a simple way much of the 
physics of the solid phases of soft core particles. Further- 
more, we believe the present theory is of general interest 
to the classical DFT community, since many DFT the- 
ories are unable to describe solid-solid coexistence. For 
example, in Refs. [2lL |22( the authors applied DFTs that 
are very successful for hard-spheres to fluids composed of 
particles interacting via Yukawa and inverse power pair 
potentials and found them not to predict the bcc phase. 
The present work may give some insight into what is 
required in a DFT in order to describe a solid-solid tran- 
sition - see also Ref . [23| . 

This paper is laid out as follows: In Sec. [H] we describe 
the DFT theory. In Sec. IIIII we apply the theory to the 
GCM and then in Sec.HVIto the YCM. Finally, in Sec. 
0we discuss our results and draw some conclusions. 

II. THE DFT THEORY 
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Given an expression for the Helmholtz free energy of a 
system, one can obtain all other thermodynamic quanti- 
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ties. It can be shown that the Helmholtz free energy T 
is a unique functional of the one body density profile of 
the system, p(r) [jj. We can divide the Helmholtz free 
energy into two parts: 

T\p\=T*\p\+TM. (3) 

The first term is the ideal gas contribution 0: 



Tid\p\ = k B T J drp(r)[ln(p(r)A 3 ) - 1], 



(4) 



where T is the temperature, k B is Boltzmann's constant 
and A is the thermal de-Broglie wavelength. The excess 
part is given formally by [|| : 

Fex[p} = ^ J dr J v(r,r')p(r)p(r') / dag(v,v';a) 

(5) 

where g(r, r'; a) is the inhomogeneous radial distribution 
function corresponding to a system of particles interact- 
ing via the pair potential av(r, r'), i.e. Eq. 10 is formally 
derived by 'turning-on' the interactions between the par- 
ticles (via the parameter < a < 1) and integrating 
gjr, r'; a) as a is increased from to 1 keeping p(r) fixed 
PI- The main approximation in our theory involves re- 
placing the function dag(r, r'; a) by a simple ansatz. 
This is given below; first we make a few remarks about 
P(r). 

When in the uniform fluid state, the one-body density 
is a constant, p(r) = p = N/V, the average number 
density, where N is the number of particles and V is 
the volume of the system. However, when the system 
freezes into a solid the density becomes periodic, i.e. the 
symmetry breaks and p(r) = p(r — Rj), where R^ is a 
lattice vector for the solid phase. An approximation that 
is often made in DFT studies of freezing 0, [24|, is to 
assume that the density profile of the solid is made up of 
Gaussian peaks centred on each lattice site: 



JY 



p(r)=^G(r-R 2 



(6) 



with 



G(r-Ri) = (^^expI-a^-R,) 2 ], (7) 

where a is a parameter which describes how localised the 
particles are around each lattice site. This density pro- 
file assumes the normalisation condition that there is one 
particle per lattice site. Of course, this assumption need 
not necessarily be true. We can expect to find vacan- 
cies in the crystal, and indeed for the present soft-core 
systems we may find another type of defect: double occu- 
pancy, since the particle cores can overlap. However, we 
expect the proportion of defects to be small, and assume 
a perfect crystal with sites singly occupied. 



Together with the assumption that p(r) takes the form 
in Eq. © we make the further assumption that T ex takes 
the form: 



J dr J dr' «(r - r')G(r - R 4 )G(r' - Rj). 

(8) 

This constitutes an RPA like approximation for 
J dag(r,r';a) in Eq. (J5J, with the 'self-energy' term in 
the summation over lattice vectors, R^ = Rj, being sub- 
tracted. When a is sufficiently large so that there is neg- 
ligible overlap between the density peaks on neighbouring 
lattice sites, our assumption is equivalent to: 



/ dasr(r,r';a) 
Jo 




(9) 



where the length I ~ bi/2 and b\ is the distance be- 
tween nearest neighbour lattice sites. This approxima- 
tion therefore constitutes quite a drastic simplification of 
the function J* dag(r, r'; a), which neglects much of the 
information about correlations in the system that this 
function contains. 

Eq. © is very appealing because it takes the form of 
a double convolution and can therefore be written in the 
form: 



1 



1 
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dkexp(ik-Ry)u(fc)G(fc)G(£;) 

(10) 

where R^ = Rj — R; and G(fc) is the Fourier transform 
of G(r). Since G(r) is a Gaussian function, then so is 
G(fc). 

The ideal gas contribution to the Helmholtz free energy 
flUl also takes a simple form if we assume p(r) is given 
by Eq. © and we further assume that a takes values 
sufficiently large that the overlap between the Gaussian 
density peaks on neighbouring lattice sites is negligible. 
Then the ideal gas part of the Helmholtz free energy, Eq. 
(gj, is simply (see e.g. plf) 



F id {T,N,a)=Nk B T 



2 \ 7T 



(11) 



Given a pair potential v(r), for which the Fourier trans- 
form v(k) exists, Eqs. (fTUfl and (fTTft together provide 
an expression for the free energy T which is a function 
of temperature, the average density p, the parameter a 
and the set of lattice vectors {R^}. For a given state 
point (T, p) and lattice structure, we assume that the 
parameter a is determined by the minimisation condi- 
tion (dJ- '/ 1 da) a=amin = 123 and we assume that the 
Helmholtz free energy F = !F(a m i n ). We can therefore 
calculate the Helmholtz free energy for a number of can- 
didate crystal structures and then the equilibrium crystal 
structure is that with the lowest free energy. In order to 
determine the melting phase boundary we could compare 
this minimal value of f — F/V with that calculated from 
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the theory applied to the liquid state and then perform 
the common tangent construction between these free en- 
ergies (which is equivalent to equating chemical poten- 
tials and pressures in the coexisting phases 0). In the 
liquid state, where p(r) = p, the ideal gas contribution 
to the Helmholtz free energy (J3J becomes 



T id {p) = Nk B T[\n(pA 3 ) - 1], 



(12) 



and the excess contribution, obtained from Eqs. (J5J) and 
@, is: 



T^{p) = 2nNp / drr 2 v(r) 



(13) 



The latter constitutes a very crude approximation for the 
liquid state free energy. Given that our aim here is to 
construct a theory which is above all simple, but which 
is still able to provide a qualitative description of the 
solid phases of soft-core particles, we choose to employ 
a Lindemann criterion |2y| to calculate the solid melt- 
ing curves, rather than compare our solid free energy 
with that of the liquid. The Lindemann criterion sim- 
ply states that when the root-mean-square displacement, 
a = ((r 2 ) — (r) 2 ) 1 / 2 , of a particle about its equilibrium 
position is roughly 10% of the nearest neighbour distance 
&i, the crystal will melt. For the Gaussian density pro- 
file, Eq. 10, a = y/3/2a and we determine approximate 
melting boundaries from the locus defined by a/b± =0.1. 

One can improve upon Eq. © as an approximation for 
the density profile in the crystal: in order to incorporate 
the effects of anisotropy in the density peaks around each 
lattice site one can assume the density profile is of the 
following form: 



N 



3/2 



= E - e^ r - R *) [1 + ™ 2 /(r - R*)], (14) 



y 4 + z A — 3r 4 /5 is the 



where the function /(r) = x A 
lead ing term for the unit cell anisotropy in cubic lattices 
and r is an anisotropy parameter. One then 
minimises the free energy with respect to r, as well as a. 
We attempted such an approach for the GCM, but we 
found it made no significant change to the phase diagram 
that we obtained from the simple choice ©. 



III. APPLICATION TO THE GCM 

In the GCM, the interparticle pair potentials are given 
by Eq. ©. The Fourier transform of this Gaussian po- 
tential is also a Gaussian and so the excess Helmholtz 
free energy given by Eq. (|T[j|) takes the particularly sim- 
ple form: 



GCM 



E 

^3 



3/2 



2A 3 



• exp(- 



-iRl), 



(15) 



R; 



where 7 = (1/A 2 + 2/a)" 1 and R. K 
(|llf> and (|15[) as our approximation for the Helmholtz free 



Using Eqs. 
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FIG. 1: Phase diagram of the GCM. The solid line is the 
bec melting boundary and the dashed line is the fee melting 
boundary, both determined using the Lindemann criterion - 
see the text. The dot-dashed line is the locus of points where 
the Helmholtz free energies of the bee and fee phases are equal. 
The dotted lines are the phase boundaries obtained from the 
computer simulations of Prestipino et al. |l9f . 



energy of the crystal, we calculate the phase diagram for 
the GCM. The results are displayed in Fig. ^ Recently, 
Prestipino et al. |19| made an accurate determination of 
the GCM phase diagram using Monte Carlo simulations, 
so we are able to compare with these essentially exact 
results (see also Ref. [l8|). We find, as did Prestipino et 
al. |19| , that at low temperatures, on increasing the den- 
sity, the fluid first freezes to form a face-centred-cubic 
(fee) crystal, and then on further increasing the density 
there is a transition from the fee to a body-centred-cubic 
(bec) crystal - see Fig. Performing the common 
tangent construction between the bec and fee free ener- 
gies, we find that the two-phase region between the two 
crystal phases is very narrow, the difference between the 
coexisting densities ApA~ 3 ~ 10~ 4 - see also the inset 
to Fig. 9 of Ref. ^^|. Since we are mostly interested in 
providing a simple theory which accounts for the topol- 
ogy of the phase diagram, we determined the density at 
which the Helmholtz free energy of the be e eq uals that 
of the fee structure for a given temperature [3l| . The re- 
sulting line is plotted in Fig. ^ The present theory also 
accounts for the most striking feature of the GCM phase 
diagram: the re-entrant melting of the bec phase - i.e. for 
a given (low) temperature, on increasing the fluid den- 
sity it freezes, but on further increasing the density the 
crystal remelts. The high density phase of the GCM is a 
fluid. This means there is also a maximum temperature 
for which there is crystal. The present theory predicts 
this maximum to be at a temperature kgT/e ~ 0.012, 
whereas it is actually at fc^T/e ~ 0.009 |19|. In general, 
the present theory over estimates the region of stability 
for the solid phases. 

We also determine an approximate melting boundary 
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FIG. 2: Same as Fig.0 except now the solid line is the locus of 
points where the Helmholtz free energies of the bcc and liquid 
phases are equal. The liquid state free energy is calculated 
using the rather crude approximation in Eq. (1131 . 



This bending over of the fcc-bcc boundary was ascribed 
to anharmonic effects in the crystal |Tflj . In an attempt 
to incorporate this effect we assumed the density took the 
form given in Eq. (|14|) . Within the present theory anhar- 
monic effects do indeed become more significant in the 
high temperature part of the fee portion of the phase di- 
agram. However, they do not become significant enough 
to result in the fcc-bcc boundary curving over, as in the 
results of Prestipino et al. 19]. Assuming l|14(l leaves 
the phase diagram unchanged when plotted on the scale 
given in Fig. ^ 

For the results presented in Fig. ^ we summed over 
40 shells of lattice vectors (including any more does not 
change the value calculated for the free energy). How- 
ever, it is interesting to note that if one includes only the 
first two shells - i.e. nearest neighbour and next nearest 
neighbour contributions only, then the results are quali- 
tatively unchanged |32|. This is not too surprising, given 
the short ranged nature of the GCM pair potential. 



IV. APPLICATION TO THE YCM 



by calculating the locus in the phase diagram where the 
bcc free energy equals the liquid state free energy [3l| . 
where the liquid Helmholtz free energy is calculated using 
the crude approximation, Eq. (|13fl . which for the GCM 
becomes: 



2X1 



cxp(-A 2 Z 2 ) +erfc(A0 



( 16 ) 

where erfc(x) = 1 — 2tt 1 / 2 J* dt exp(— t 2 ) is the com- 
plimentary error function. Whilst we know I ~ b\/2, 
there is no constraint on a particular value. We choose 
the value I ~ 0.58(2/p) 1 / 3 (recall that for the bcc crys- 
tal bi = {VS/2)(2/p) 1 ^ 3 ). This value of I is chosen so 
that the maximum temperature that freezing occurs is 
predicted to be roughly at the same temperature as that 
from the results of Prestipino et al. 19]. Note, on the 
scale of Fig. O using this approach there is no differ- 
ence between the predicted locations of the bcc-liquid 
melting boundary and the fcc-liquid melting boundary. 
Given the crude nature of the theory for the liquid, the 
results are surprisingly good - see Fig. 

The main feature of the GCM phase diagram that is 
not accounted for by the present theory, and which can 
be seen in the results of Prestipino et al. ( see a l so 
Fig- DJ, is that the transition line between the fee and 
bcc phases is not the (almost) straight line predicted by 
the present theory. Instead, on following the fcc-bcc 
boundary as the temperature is increased, the bound- 
ary obtained in simulation curves over to lower densities 
and in fact reaches a maximum and then bends down 
to lower temperatures and densities, such that there is a 
small window in the phase diagram where for tempera- 
tures around fceT/e ~ 0.0035, on increasing the density, 
the fluid first freezes to form the bcc, then the fee and 
then the bcc phase again, before finally re-melting [l9j . 



The Fourier transform of the YCM pair potential, Eq. 
(JTJ, is v(k) = 47reA _1 (A 2 + k 2 )^ 1 . Using this, together 
with Eq. IjlQI) , we obtain the following expression for the 
YCM excess Helmholtz free energy: 



T. 



YCM 



\ - ec 



A 2 /2q 



4Ai? 4 



e -XR Werfc 



-e A ' Ri3 erfc 



A 



2a 



(17) 



This approximation, together with Eq. H 1 1 p . is our ex- 
pression for the Helmholtz free energy of the solid phases 
of the YCM. For a given state point (p, T) and set of 
lattice vectors {Ri}, we minimise the free energy with 
respect to the parameter a. We estimate the melting 
boundaries using the Lindemann criterion - i.e. the locus 
defined by a jb\ = 0.1. In Fig. we display the resulting 
phase diagram. For the YCM, the crude estimate for the 
liquid state excess Helmholtz free energy, Eq. I|13|) . does 
not give physically acceptable results. For some choices 
of I in (|13f) it (incorrectly) predicts that the YCM ex- 
hibits re-entrant melting. We therefore employ only the 
Lindemann criterion for determining melting boundaries 
in this section. 

The low density portion of the YCM phase diagram is 
qualitatively similar to that of the GCM - i.e. for suffi- 
ciently low temperatures, on increasing the density the 
fluid first freezes to form an fee crystal, then, at higher 
densities, there is a transition to the bcc. The biggest dif- 
ference between the YCM and the GCM phase diagram 
is that there is no re-entrant melting in the YCM. This 
is because the divergence of the YCM pair potential as 
r — > means that the particles behave more and more 
like hard spheres as the density is increased, where the 
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FIG. 3: Phase diagram of the YCM. The solid line is the 
bec melting boundary and the dashed line is the fee melting 
boundary, both determined using the Lindemann criterion - 
see the text. The dot-dashed line is the locus of points where 
the Helmholtz free energies of the bec and fee phases are equal. 
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FIG. 4: Same as Fig.|3 except here the results are plotted in 
terms of the variables k and F. The symbols joined by solid 
lines are the simulation results of Hamaguchi et al. |34| |. 



effective hard sphere diameter depends upon the state 
point (T,p), so that there is freezing of the YCM at all 
temperatures. For the GCM this is not the case. Note 
that a divergence in the pair potential at r — > does not 
automatically mean that there is no re-entrant melting 

m 

Just as for the GCM, the present theory for the YCM 
also fails to account for the curving over to lower densities 
of the fcc-bcc boundary as the temperature is increased 
in the YCM. This can be seen in Fig.0| where we plot the 
YCM phase diagram that we obtain together with that 
of Hamaguchi et al. (see also Refs. jM H El S 
I39l l40ll41j ). plotted in terms of the variables T = f3e/\a 



and k — Xa, where a = (3/47TP) 1 / 3 . These are commonly 
considered variables when using the Yukawa potential 
to model the interactions in plasma systems. We see that 
the present simple theory is able to account qualitatively 
for the YCM phase diagram. 



V. DISCUSSION AND CONCLUSIONS 

We have constructed a simple DFT for the solid phases 
of soft-core particles. The theory is able to account for 
the transition from a fee to a bec solid and also able to 
determine whether the system exhibits re-entrant melt- 
ing, as is the case for the GCM, or not, as in the case 
for the YCM (when a Lindemann criterion is used to de- 
termine the melting boundaries). This makes the DFT 
useful, since many DFTs are not able to describe solid- 
solid coexistence in soft core fluids 0, ^ . There are 
some DFT theories able to describe solid-solid transi- 
tions in other (hard-core) model fluids - see for example 
Refs. I n fact, the present theory bears some 

similarities in its structure to that of Likos et al. 0, ^2 • 
This can be traced to the use of the Gibbs-Bogoliubov 
inequality to construct their theory. When one applies 
this inequality, one obtains the following equation (Eq. 
(4.18) in Ref. Q): 

F\p\ < MP] + \ J dr / dr' 5o (r,r>(r)p(r>(r - r'), 

where go and T§ are the pair distribution function and 
the Helmholtz free energy respectively of the reference 
system; Likos et al. 0, (43 used a hard-sphere fluid as 
the reference system. If we compare Eq. (I18|) with Eq. 

we can see that depending on our approximations for 
<7o(r, r') in (|T%)l and L dag(r,r';a) in ||SJ, one can end 
up with theories that have a similar structure. 

One also sees similar features when we compare our 
expression for the GCM Helmholtz free energy, Eqs. J ll |) 
and (|15|) , with that obtained by Lang et al. [l], [l8| us- 
ing the Gibbs-Bogoliubov inequality together with the 
Einstein model as their reference system. Their result- 
ing Helmholtz free energy is almost exactly the same as 
that in the present theory (the free energies differ only by 
ksT per particle, with that of Lang et al. being lower). 
In the theory of Lang et al. Q, Il8| the Einstein model 
spring constant is the variational parameter for minimis- 
ing the free energy, whereas in the present theory it is 
the parameter a in Eq. ©. However, formally these 
parameters play exactly the same role. Use of the Gibbs- 
Bogoliubov inequality therefore seems to lead to theories 
with a structure similar to the theory presented here. 
When taking our present approach, i.e. using the Linde- 
mann criterion to determine the melting boundaries, the 
difference of NksT between the Helmholtz free energy 
of the present theory and that of Lang et al. 0] makes 
no difference since this term is independent of a. How- 
ever, it would matter if we were to compare our result 
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for the solid free energy with that obtained for the liquid 
from some other theory, more accurate than that from 
Eq. P|. 

Some soft core fluids exhibit freezing to states with 
multiple occupancies of each lattice site 0, E3- I 11 or ~ 
der to apply the present theory to such systems, some 
modification of the theory is required. Firstly, one must 
assume a generalisation of Eq. JBJ for the density profile 
of the crystal: 



N 



3/2 



-a(r-Ri) 



(19) 



where rj is the average lattice site occupancy. 77 should 
be treated as a parameter to minimise the Hclmholtz free 
energy, in the same way as with the parameter a. In this 
case, there would be two minimisation conditions to be 
satisfied: {dT / da) a=amin = and {d!F / dn) v=rimin = 0. 
One would then assume that the Helmholtz free energy 
F = T{a m in,ri rn i n ). This would also be the scheme to 
apply if one intended to study the effect of lattice defects 
in the present YCM and GCM systems. However, in 
these cases one would expect r\ ~ 1 . For systems exhibit- 
ing multiple occupancies of each lattice site, we expect 
one would also have to make a different approximation 
for the function J* dag(r, r'; a) in Eq. JSJ. We would 
propose the following generalisation of Eq. JSJ: 



From a more general point of view, the present theory 
seems to provide a good qualitative description of the 
solid phases of soft-core systems. The Yukawa potential, 
(Q, is used to model the effective interaction between 
charged colloidal particles 0, B> 01 G3 • For example, the 
phase diagram of polystyrene particles suspended in a 
potassium chloride solution can be mapped on to that of 
the YCM 0| . The present theory should also be relevant 
to soft matter systems, for example polymeric micelles 
|22|. star polymers [33| and dendrimers |46|. Given an 
effective pair potential v(r) between such objects pj, one 
could use the present theory to calculate an approximate 
phase diagram. 
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